COMPUTATION OF ENGINE NOISE PROPAGATION AND 
SCATTERING OFF AN AIRCRAFT 


J. Xu*, D. Stanescu*, M.Y. Hussaini*, and F. Farassat^ 

*School of Computational Science and Information Technology, Florida State University 

Tallahassee, Florida 32306-4120 

t NASA Langley Research Center, Hampton, Virginia 23681-0001 


Abstract 

The paper presents a comparison of experimental 
noise data measured in flight on a two-engine busi- 
ness jet aircraft with Kulite microphones placed on 
the suction surface of the wing with computational 
results. Both a time-domain discontinuous Galerkin 
spectral method and a frequency- domain spectral el- 
ement method are used to simulate the radiation of 
the dominant spinning mode from the engine and 
its reflection and scattering by the fuselage and the 
wing. Both methods are implemented in computer 
codes that use the distributed memory model to 
make use of large parallel architectures. The results 
show that trends of the noise field are well predicted 
by both methods. 

INTRODUCTION 

The fan inlet and exhaust noise represents one 
of the major components of the noise signa- 
ture of an aircraft at take-off and landing. The 
noise radiated to the far field by the engine 
of an aircraft is largely influenced both by the 
flow around the wing and fuselage and by the 
scattering from various other aircraft surfaces. 
In principle, it is possible to reduce the air- 
craft noise footprint by taking advantage of en- 
gine and wing location and manipulating the 
flow around the aircraft. Experimental inves- 
tigations of these phenomena are difficult to 
perform and extremely expensive. Numerical 
simulations offer a relatively inexpensive alter- 
native, and such simulations are becoming in- 
creasingly attractive due to the recent advances 


in both computer architecture and computa- 
tional methods. To date, most measurements 
and modeling of engine noise are confined to 
isolated engines [1, 2, 3]. However, Stanescu, 
Hussaini, and Farassat [4] have recently com- 
puted the engine noise propagation and scat- 
tering for a generic aircraft configuration by nu- 
merically solving the Euler equations by a dis- 
continuous Galerkin spectral element method. 
The recent popularity of such methods in aero- 
dynamic applications stems from the fact that 
they require relatively fewer points per wave- 
length, they have low dispersion and dissipation 
errors [5, 6, 7], they have geometric flexibility, 
and they are compact, robust and inherently 
parallelizable [4, 8]. In a more recent work [9], 
the same authors presented a simultaneous re- 
search initiative consisting in the development 
of a spectral element method for the solution, 
in the frequency domain, of the acoustic poten- 
tial equation in the presence of a non-uniform 
flow field. The two methods are believed to be 
complementary tools useful for prediction of the 
tonal sound field in the near field of the engine. 

The purpose of the present effort is to in- 
vestigate the feasibility of large-scale aircraft 
noise simulations, and validate them with the 
available experimental data. To that end, we 
employ the aforementioned numerical method- 
ology in the investigation of engine noise prop- 
agation and scattering off an actual two-engine 
jet aircraft. After a brief discussion of the two 
methods in the next section, we present the nu- 
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merical results obtained for several flight condi- 
tions and compare them with the available ex- 
perimental data. Although the modal composi- 
tion of the source could not be obtained exper- 
imentally and the actual flight conditions could 
not be replicated exactly in the computations 
presented here, results show that the noise sig- 
nature obtained from computations matches the 
trend of the experimental data. 


PROBLEM FORMULATION AND 
SOLUTION TECHNIQUE 

Computational model 

We assume that the engine noise source is 
known and consider noise propagation and scat- 
tering in the left half space implicitly assuming 
symmetry of the problem about the y = 0 plane 
which is chosen to bisect the aircraft along the 
fuselage. For the computational purpose, the 
half space is truncated into a computational do- 
main comprised of a bounded physical domain 
with a damping layer surrounding it. The lat- 
ter is used to ensure the physical domain re- 
mains uncontaminated by reflections. The sur- 
face that separates the computational domain 
from the surrounding medium is denoted by 
Too • This governing equations are solved in non- 
dimensional form. The reference quantities for 
non-dimensionalization are: poo for th® density, 
Coo for the velocity components, PooC^ for the 
pressure, the radius R of the noise-source disk 
for distance, and — for time. The total domain 
of computation (including the damping layers) 
is defined in the non-dimensional Cartesian co- 
ordinates as 21.4 < X < 40.0, —12.0 < y < 0.0 
and —1.8 < 2 : < 11.2. The computational do- 
main with the embedded aircraft is depicted in 
Fig. 2. As the propagation distance is rela- 
tively small, viscous effects are neglected and 
the problem is assumed to be governed by in- 
viscid compressible flow equations. 

The computational domain is covered by a 
grid of non-overlapping general hexahedral ele- 
ments that can have curved boundaries. The 
ICEMCFD Hexa commercial package is used 
to generate the unstructured hexahedral grid 


around the aircraft configuration. Once an un- 
structured grid of hexahedra is generated, an 
attractive new feature of this package allows for 
the generation of points along each of the edges 
of the hexahedral mesh, which can be either a 
Legendre-Lobatto or a Chebyshev-Lobatto dis- 
tribution for a specified polynomial of degree N . 
Fig. 3 shows the hexahedral representation of 
the underlying geometry with Gauss-Legendre 
point distribution, where N = 5. All the neces- 
sary point coordinates can then be computed by 
interpolation based on the spectral interpolants 
along the edges (to obtain coordinates at the 
Gauss points from the Lobatto points on the 
edges) followed by three-dimensional transfinite 
interpolation [10] on the faces and inside the el- 
ements. 


Boundary conditions 


A zero normal velocity boundary condition is 
imposed on the symmetry surface y = 0 (this 
supposes that the engines are symmetrically 
placed on either side of the fuselage and rotate 
in opposite sense) as well as on the fuselage, 
nacelle and wing surfaces. The boundary condi- 
tions on the other sides of the computational do- 
main that make up Too are treated by a damp- 
ing layer method [11]. The damping layer is 
about 3.5 wide in the x—,y— and 2 — directions. 
Waves incident on this layer are damped and 
reflections into the physical domain of interest 
are minimized. This is obtained by modifying 
the governing equations through the addition of 
a damping term in the form 


dQ 

dt 


+ V-T 


-cr(x)Q 


( 1 ) 


where the damping parameter is made to vary 
from 0 at the interior limit of the damping layer 
to a maximum value on Foo according to a power 
law 


<t(x) = (TmY^ 

i 


( X,- 

* - a:*"* 


/3 


( 2 ) 


with a;™* and xf^^ the coordinates of the interior 
and exterior limits of the absorbing layer, limits 
that lie along planes on which one coordinate is 
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constant. 

The engine tone noise source is specified as 
a combination of spinning modes on a surface 
r f conveniently situated inside the nacelle. For 
a single spinning mode with azimuthal order s 
and radial order d, usually denoted as (s,d), 
the perturbation of the flow variables from the 
mean flow quantities (denoted by bars) is given 
by [11, 12]: 


/ p-p \ 

P-P 
Vx - Vx 

Vr — Vj. 

\ ve-ve J 


/ 
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Emikmdr) cos G ^ 
^ • Ern{kmdr) cos O 
COS 0 

t:$-E'rn{kmdr) Sine 
j^^-Em{ksd'r)cose ^ 


where 0 = kxX + mO — u)rt, kmd = 

{uirlcY — kx'^ , and uir = The function 

Emik'/nd'c) — Jmikmd'C^ E (lYrnikrndE) is the duct 
eigenfunction with Jm and Ym denoting the 
Bessel functions of the first and second kind, 
respectively. The noise source is specified on 
the circular disk centered at (34.7,4.6,5.3). 


Time domain formulation 


For the time domain formulation, the governing 
equations are considered the Euler equations in 
Cartesian conservation form, 


dQ ^dFd 
dt ^^dxd 


0. 


(4) 


where the state vector Q and the flux vector Ed 
are given by 


representing the solution in each element by 
spectral basis functions defined on the interval 
[—1,1]. Under the mapping, Eq. (4) becomes 


dQ ^dFd 
dt ^^d^d 


0 , 


(6) 


where Q and F are the transformed components 
of the state and flux vectors 


Q = JQ, Fd = Jj2p^F^, (7) 

m=l"Xm 

and J is the Jacobian of the transformation. 
The computational space coordinates are de- 
noted by either (<ti, (^ 2 > <? 3 ) or (6h>C) hereafter 
for convenience. 

Let the space of polynomials of degree N in 
^ G [—1, 1] be denoted by Pn- A basis for this 
space can be constructed using the Lagrange 
interpolating polynomials hj,j = 0,1, . . . , N, 
through the A1 + 1 Gauss-Legendre [13] quadra- 
ture nodes i = 0,1, . . . , N . A discontinuous 
Galerkin approximation is obtained by requir- 
ing 

{Qt, 4>ijk) + (Vg • F, = 0 (8) 

where (•, •) represents the usual inner prod- 
uct, and 4>ijk = hi{^)hj{r])hk{C) are the basis 
functions of P^. 

Using the divergence theorem and Gauss 
quadrature, expanding the boundary integral 
and performing some algebraic manipulation, 
the final discrete form of the equations govern- 
ing each variable at the Legendre-Gauss points 
are given by 


^ p \ 


( PVd \ 

pvi 


pviVd +p5id 

pV2 

, Fd = 

pV2Vd + p52d 

pV3 


pvsVd + phd 

\ pE J 
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with p the fluid density, E the internal energy, 
p the pressure, and Vd{d = 1,2,3) the velocity 
components. 

Each element of the grid is mapped onto the 
master element Q,m = [— 1, 1]^ with an isopara- 
metric transformation for the expediency of 


F, (9) 

where the right-hand side is a sum of discrete 
differential operators acting on the flux values 
of an element, which include values on the ele- 
ment faces. Here, the differential operator D^F, 
for example, is defined as 
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where F* denotes a common face flux, which 
can be computed directly from the readily avail- 
able values of the state vector. D^F and D^F 
follow by obvious permutations. 

As the solution is approximated by a poly- 
nomial that passes through interpolation nodes 
distributed within the elements, a mismatch en- 
sues when the interpolants are evaluated at el- 
ement interfaces. This mismatch in the solu- 
tion at element boundaries is resolved by solv- 
ing the Riemann problem for the flux there (just 
as in the flnite volume method) [14, 15]. This 
leads to a semidiscrete form of Euler equations, 
which is simply an ordinary differential equa- 
tion (ODE) system. The resulting ODE sys- 
tem is integrated in time using a low-storage 
Runge-Kutta scheme optimized for wave prop- 
agation [16]. Acoustic perturbations are ob- 
tained at each time step by subtracting the 
mean flow from the total flow variables, and 
the RMS pressure is obtained by integrating in 
time the acoustic pressure. This integration is 
only performed after sufficient time is allowed 
for the acoustic signal from the source to prop- 
agate through the computational domain and 
establish a periodic acoustic ffeld. 

Frequency domain formulation 

The equation governing the acoustic field is in 
this case obtained by considering the flow irrota- 
tional, so that the continuity equation becomes 

|j+V.(pV$) = 0 (11) 

where p is the fluid density, and d> is the to- 
tal velocity potential, related to the velocity by 
V = Vd*. Under the isentropic assumption, the 
momentum equation is reduced to an algebraic 
relation relating the density to the velocity po- 
tential as 


1 


1 - (7 - 1) 



(v$)^-M^y 


7-1 


( 12 ) 


where M^o is the far field Mach number and 
7 the specific heats ratio. Consider the un- 


steady flow field resulting from the superposi- 
tion of small acoustic perturbations, denoted by 
a prime, on a steady mean flow denoted by an 
overbar: p = p + p' and T = T + T'. The par- 
tial differential equation governing the acoustic 
perturbations is 

^ + V • (pV$' + p'Vi) = 0 (13) 

with the following relation relating the acoustic 
density to the acoustic velocity potential, ob- 
tained by linearization of equation ( 12 ): 


/ 

P 







(14) 


For a frequency domain approach, the acoustic 
potential is considered to be of the form = 
4>(x,y,z)ex]){iujt). In view of a weak formula- 
tion, the governing equation (13) is multiplied 
by a test function T' = y, z) exp{—iu>t) and 
integrated using the divergence theorem to yield 

[ i>^rcin - / VV- • + p'W^)rdQ+ 

Jq ot jQ 

I il}{pS/(f) + p'^$) ■ nrds = 0 

(15) 

On the aircraft surface both the mean flow and 
the acoustic normal velocity component are sup- 
posed to be zero, so the surface integrand can- 
cels there. Furthermore, the integrand is set to 
zero artificially on the Too surface, which does 
not sensibly affect the computed solution since 
the acoustic field is strongly damped in the ab- 
sorbing layer. 


Let Z denote the complex vector space of 
functions that are continuous on 0 , whose re- 
strictions to an element are polynomials of de- 
gree at most N in each variable, where N is 
a specified integer, and Zpj. C Z the subset 
of functions that satisfy the Dirichlet boundary 
condition on the source surface F j. Substituting 
the expression of the density from the linearized 
momentum equation, the following variational 
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problem is obtained: find (j) £ '^D such that 


trix A can then be written in the form 


/ ^ + iuju + 

Jnc L 
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(U^ - C?) (j)y-tljy + 

{up - C^) 4>zpz + UV {(pxi’y + 4>yi’x) + 

UW {4>xPz + Pz'Px) + VW {(f^ylpz + 4>z'4’y) ~ 

pa {iujp + upx + vpy + w<l>z)] dPt = 0 

(16) 

holds for any p G Zpj. . 


The previous equation is discretized by a 
Chebyshev spectral element method. To this 
end, a basis for Z is constructed using tensor 
products of the Lagrange interpolants through 
the Chebyshev-Gauss-Lobatto points in the ele- 
ment. Upon evaluating the integrals, a complex 
linear algebra problem of the form A{p} = {6} 
is obtained, where {p} is the vector of point 
values of the complex-valued acoustic potential 
p. The solution of this system is obtained us- 
ing a Schur complement domain decomposition 
method implemented using the Message Pass- 
ing Interface (MPI) standard. The matrix is 
stored in sparse mode (i.e. only the non-zeroes 
are stored), with each processor only storing a 
number of lines in the matrix. Let us denote by 
P the total number of processors, the computa- 
tional domain 0 is subdivided in as many parti- 
tions, and the unknowns situated on the surface 
B which separates the partitions are numbered 
last in the system. For every processor p, there 
will be a number of unknown p values located 
on B. The vector of unknowns is partitioned as 

{P] = 4>b] (17) 

where p^j denotes all the unknowns in subpar- 
tition p not located on B. The right-hand side 
vector {6} is partitioned accordingly. The ma- 


and straightforward elimination of the terms be- 
low the main diagonal leads to 



(19) 

whereas = fcs— The problem 


p 

has thus been reduced to solving a reduced sys- 
tem with matrix S = Abb — 

p 

for the points on B only, followed by a solution 
on each domain of the interior problem. The 
matrix S is much denser than the original ma- 
trix A and its direct computation and storage is 
not efficient or even possible. However, for an 
iterative method, only the action of S' on a vec- 
tor is needed, and once the sparse, distributed, 
matrix Abb formed, this action can be com- 
puted by matrix-vector multiplications and so- 
lutions with A^j which are local operations on 
processor p and do not require communications, 
followed by accumulation in the global vector 
pB- All computations can be conveniently im- 
plemented by use of the high level primitives in 
the PETSc [17] package for efficient solution of 
partial differential equations. 

RESULTS AND DISCUSSION 
Experimental data 

The data available for comparison was collected 
in Moo = 0.3 flight at 500 foot above ground 
level. The average Mach number at the source 
disk based on the mass flux through the engine 
is approximately Mj = 0.53. The Blade Pass- 
ing Frequency (BPF) tone was measured us- 
ing several Kulite microphones located on the 
wing suction side at different angles from the 
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nacelle axis as shown in Fig. 1, at a relatively 
high power setting of the engine. Modal com- 
position of the signal could not be satisfactorily 
established, however there were indications that 
the dominating mode is the spinning mode with 
azimuthal order m = 22, and m= 18 is also ex- 
pected to have a large contribution. The data 
is presented as Sound Pressure Level (SPL) val- 
ues, normalized such that the value at 20 de- 
grees corresponds to zero SPL. 

Time domain results 

Propagation of both modes was simulated in 
the time domain with a quiescent tnedium sur- 
rounding the aircraft in a first attempt to de- 
termine the dominating mode and conduct a 
first test of the method on this geometry. The 
BPF tone corresponds to a reduced frequency 
ujr = 26.3. Propagation and radiation of modes 
(18,0) and (22,0) was computed separately us- 
ing an unstructured grid with 103, 105 elements. 
The solution is approximated by a sixth-order 
Legendre polynomial in each element, raising 
the number of Gauss-Legendre discretization 
points to 22,270,680 in the computational do- 
main that includes the damping layers. The 
computations used one node (32 processors run- 
ning at l.lGHz) of an IBM Regatta-type SP4 
machine and each one lasted about 10 days. An 
arbitrary value has been used for the ampli- 
tude of the incoming mode, which is not known 
from the experiments. Therefore, for the pur- 
pose of data comparison, computational data 
was matched with experimental data at the 60 
degree microphone location, where the peak in 
SPL was noticed experimentally. 

Fig. 4 shows a snapshot of the acoustic pres- 
sure contours on the surface of the aircraft at 
non-dimensional time t = 44, immediately be- 
fore starting integration for the RMS pressure 
computation, for spinning mode (18,0). The 
computation indicates that radiation of mode 
(22, 0) has a pattern that is completely differ- 
ent from the experimental data, presented in 
Fig. 5. However, mode (18,0) seems to pro- 
duce a SPL distribution on the wing with the 
same characteristics as the experiments, Fig. 6. 


The quantitative difference in levels may stem 
from the mean flow effect. Indeed, increasing 
the mean flow Mach number in the duct from 
its zero value used for the present results will 
determine the mode to be more cut-on, with an 
immediate effect that the main radiation lobe 
will hit the wing at a lower angle location. Thus, 
SPL levels at angles lower than 60 degrees are 
expected to increase, while decreasing at 70 de- 
grees. The effect of the mean flow, on the other 
hand, is expected to be in the opposite sense, 
but not as strong. It must also be considered 
that other sources of noise, not modeled in this 
computation (i.e. airframe noise) will increase 
the value of the measured SPL data, and this 
effect will become more important at larger dis- 
tances from the engine (lower angles in the fig- 
ures) . 

Frequency domain results 

The effect of the mean flow has not yet been 
attested with the time domain formulation, as 
our efforts concentrated lately on the develop- 
ment of the frequency domain code. The latter 
incorporates a solver for the mean flow around 
the given configuration (at this time based on an 
incompressibility assumption to avoid nonlinear 
iteration). The frequency domain code however 
has far larger memory requirements. To date, 
our largest computation with this code was done 
on the same mesh with 103, 105 elements but 
with quartic elements. There are thus a total 
of 6, 746, 736 discretization points in the mesh, 
which would generate a complex matrix with a 
total of 1.4 billion non-zeroes whose only stor- 
age in sparse mode necessitates approximately 
20GB memory. This resolution does not allow 
us to handle the exact experimental conditions, 
in particular at larger distances from the na- 
celle where the mesh becomes coarser and the 
necessary number of points per wavelength is 
not reached. Therefore, we considered a test 
case in which the mean flow around the aircraft 
was modeled based on a far field Mach num- 
ber Moo = 0.1 and a fan face Mach number 
Mf = 0.2. The problem was run on 192 pro- 
cessors of the same IBM machine, and conver- 
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gence more than three orders of magnitude was 
achieved after 2.5 days. This clearly demon- 
strates the need for a better multi-level precon- 
ditioner for the reduced matrix S. We must 
point out, however, that to our knowledge this is 
the largest problem of this type reported in the 
literature. Moreover, solving the complex linear 
algebra problem associated with discretization 
of Helmholtz-type equations is a well-known is- 
sue in the numerical analysis and parallel com- 
puting community and no satisfactory solution 
has, to our knowledge, been found as of now. 
Results obtained from this simulation for mode 
(18, 0) are presented in Fig. 7, which shows that 
the presence of the mean flow has the antici- 
pated effects as discussed above. 
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Fig. 1: Location of microphones on the wing and corresponding experimental data. 
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Fig. 2: Computational domain. 



Fig. 3: Mesh on the aircraft surface for N = 5 (quintic elements). 
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Fig. 4: Acoustic pressure contours on the aircraft surface. Mode (18,0) radiated at LOr = 26.3. 


SPL variation on the wing suction side 



Fig. 5: SPL levels on the wing surface for mode (22,0) radiated at ujj. = 26.3 in a quiescent medium. 
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SPL variation on the wing suction side 



Fig. 6: SPL levels on the wing surface for mode (18,0) radiated at = 26.3 in a quiescent medium. 


SPL variation on the wing suction side 



Fig. 7: SPL levels on the wing surface for mode (18,0) radiated at cUr = 26.3 in the presence of mean flow. 
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